arXiv:1503.00988v2 [cond-mat.soft] 27 Apr 2015 


Mean-field microrheology of a very soft colloidal suspension: inertia 

induces shear-thickening 

Vincent Demery 

Department of Physics, University of Massachusetts, Amherst, MA 01003, USA 
Laboratoire de Physico-Chimie Theorique, UMR CNRS Gulliver 7083, ESPCI, Paris, France 


April 28, 2015 


Abstract 

Colloidal suspensions have a rich rheology and can exhibit shear-thinning as well as shear-thickening. 
Numerical simulations recently suggested that shear-thickening may be attributed to the inertia of the 
colloids, besides the hydrodynamic interactions between them. Here, we consider the ideal limit of a dense 
bath of soft colloids following an underdamped Langevin dynamics. We use a mean-field equation for the 
colloidal density to get an analytical expression of the drag force felt by a probe pulled at constant velocity 
through the suspension. Our results show that inertia can indeed induce shear-thickening by allowing 
density waves to propagate through the suspension. 


1 Introduction 

Suspensions of colloids or droplets have a rich rheology [1, 2], Whereas they can have a yield stress at 
volume fractions above the glass or jamming transitions [3], they behave as simple Newtonian fluids at 
moderate densities and small shear rates. Upon increasing the shear rate, their viscosity can then decrease 
(shear-thinning) or increase (shear-thickening) [4, 5]. 

Experimentally, the suspension rheology can be investigated using macrorheology or microrheology. In 
macrorheology, a global shear rate 7 is applied, and the resulting shear stress r is measured ; the viscosity is 
then defined as 77 = r/'y [6, 7], In microrheology, the motion of a small probe in the medium is observed [8, 9]. 
Notably, in active microrheology, the probe is placed in an optical or magnetic trap and pulled at constant 
velocity v (or at constant force F) through the medium [10, 11, 12, 13, 14] (see Fig. 1). Measuring the 
drag force F on the probe (or its average velocity v ) and the corresponding drag coefficient A = F/v, 
one can use the Stokes formula to deduce the viscosity: A = 67T7?a p , where a p is the radius of the probe. 
These two approaches have their theoretical counterpart both in numerical simulations ([3, 15, 16, 17] for 
macrorheology, [18, 19] for microrheology) and analytical computations ([6, 20] for macrorheology, [21, 22] 
for microrheology). 

Shear-thinning, which is ubiquitous in experiments, is also present in most of the analytical computa¬ 
tions both for dilute [21, 22] and dense suspensions [20, 23, 24, 25, 26, 27] (see [14] for a review). It 
is commonly associated with the disruption of the equilibrium microscopic structure, which gives the solu¬ 
tion a large viscosity. On the other hand, various forms of shear-thickening exist and they are difficult to 
describe theoretically [5, 28]. While discontinuous shear-thickening may arise due to a dynamic jamming 
transition [29, 30, 31], a softer, continuous shear-thickening is induced by the formation of hydro-clusters 
due to the lubrication forces, which hold the particles together [32], However, this mechanism has a neg¬ 
ligible effect on soft-particles [33, 22, 5]: grafting polymer brushes to hard colloids can considerably delay 
shear-thickening [34], 
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Figure 1: (color online) Illustration of the system studied: a probe (red) is pulled at constant velocity through 
a dense bath of soft colloids with inertia. 

A recent numerical work addressed the role of colloids inertia on the suspension rheology, neglecting 
the hydrodynamic interactions, and showed that it can induce shear-thickening as long as the system is 
sufficiently far from jamming [15]. In this article, we provide an analytical derivation of this effect in the 
limit of very dense and soft colloids. We use microrheology to investigate the properties of this medium: we 
compute the drag force felt by a probe pulled at constant velocity through the suspension. We obtain an 
analytical expression for the drag coefficient A, which displays shear-thickening induced by inertia. 

This article is organized as follows. The model of underdamped Langevin colloids is presented in Sec. 2. 
A linearized equation for the coarse-grained density field of the bath is obtained in Sec. 3. This linearized 
equation is used to compute the response of the bath to the passage of the probe at constant velocity in 
Sec. 4. The stationnary bath density around the probe is computed (Subsec. 4.1), from which the drag felt 
by the probe is deduced (Subsec. 4.2). The analytical expressions are computed numerically and discussed 
in Subsec. 4.3. Our results are discussed in terms of the different timescales involved in the problem in 
Subsec. 4.4. Finally, the effect of the mode of driving is mentioned in Subsec. 4.5. 

2 Model 

We consider N colloids in a d-dimensional bath with positions x,(t) interacting via the pair potential V(x) 
and evolving according to Langevin dynamics, that we write in the form 

mx,(t) + A s x/(t) = Fj(t) + ri/(t). (1) 

where m is the mass of the colloids and A s the friction coefficient with the solvent. We can thus define the 
damping time associated to the inertia of the colloids as 


fi(t) is the force on the particle / and rjj(t) is the Gaussian white noise on the particle /'. The noise is 
completely defined by its correlation function 

<tj,( t)T/,(t') T } - 2T\ s 5ij5(t - 01. (3) 
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where T is the temperature (the Boltzmann constant is set to k& = 1). The force is given by the gradient 
of the potential created by the other colloids, 

Fi(t) = -Y,ViV(xi(t)- Xj (t)). (4) 

J 

The volume V of the box containing the colloids is taken to infinity, keeping the density po = A//V constant. 

A probe, which interacts with the bath colloids with the potential V p (x), is pulled at constant velocity i/ p 
through the suspension. The mode of driving, constant velocity instead of constant force, is chosen because 
it leads to easier analytical computations. The average velocity induced by a constant force applied on the 
probe has been computed without inertia for the bath particles in [27], Altough it has been shown that in 
dilute systems without hydrodynamic interactions the viscosity measured by imposing the velocity is twice the 
one measured by imposing the force [35], there is no general relation between the results obtained through 
these two modes of driving [14]. 

Without loss of generality, we can set the size of the bath colloids a, the thermal energy T and the 
friction coefficient A s to 1 (we keep the temperature explicitly in our expressions to make the role of the 
temperature clear). The dimensional counterpart of the dimensionless quantities computed below is obtained 
by multiplicating them by the appropriate factor of the size of the probe a, the thermal energy k^T and the 
friction coefficient of a colloid in the solvent A s . 

3 Linearized equation for the bath density 

Given the colloids positions x/(t) and velocities i/,(t), the bath density p(x, t ) and current j(x, t ) are defined 
by 


P(*. t) = J]5(x-x/(t)), 

/ 

(5) 

]{x, t) = J3v/5(x-x/(t)). 

(6) 


I 


For colloids with Langevin dynamics (1), Nakamura and Yoshimori derived the exact equations satisfied by 
these two fields [36]: 


d t p — -V j, (7) 

dtJ = -w - ]j)pV(\/*p) - V • J + y ^V, (8) 

t r \ p J r 

where * is the convolution product, 1 / * p(x, t ) = jV(x — x / )p(x / , t)dx' , and r/(x, t) is a Gaussian white 
noise with correlation function 


(? 7 (x, t)rj{x ', t') T ) = 2T<5(x - x')5{t - t')l- (9) 

Eqs. 7 and 8 are non-linear and contain multiplicative noise, which makes their analytical treatment 
difficult. This situation is inextricable for the real fields, but the equations can be linearized if one works with 
the coarse-grained fields p(x, t) and j(x, t), which obey [37] 


d t p = -V j, 
d t j = - 1 r 


j + TVp + pV(V * p) - yTpr] 


- V- 


JJ 

P 


( 10 ) 

( 11 ) 


T 
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In Eq. 11, we used the random phase approximation [38] to write the direct correlation function c(r) with the 
interparticle potential, c(r) = —V(r)/T. The random phase approximation is justified if the pair potential is 
weak, which is the case that we consider. 

Following [27], we can now linearize Eqs. 10 and 11 around a large homogeneous density po : writing 


p(x, t) = Po + Po / 2 0(x, t), 

(12) 

j(x, t ) = pI /2 iIj(x, t ) 

(13) 

and taking the limit po — > oo with poV(x) —> V(x), we get 

d t 4> = -V • ip, 

(14) 

d t ip = — — ip + TV0 + V(V * 0) - ^7 17 . 
r L 

(15) 

A closed second order equation can be obtained for the density field, 

rd 2 (p + d t (p = TV 2 (p + V 2 (V *</>) + vTv • r]. 

(16) 


The contribution of inertia in this expression is remarkably simple: it enters only in the first term on the l.h.s. 
Without inertia, this equation reduces to the linearized Dean equation [27], 


4 Application to a probe pulled at constant velocity 


4.1 Bath density around the probe 

In order to assess the rheological properties of the suspension, instead of shearing the material globally as in 
[15], we pull a probe particle at constant velocity v p through the medium, as in [39]. 

The interaction between the probe and a particle of the bath is given by the potential V'p(x) = V p (x)/po. 
The effect of the probe on the density field can be incorporated in the linearized equation 16 as in [27], 

rd 2 <p(x, t ) + d t (p(x, t) = TV 2 0(x, t ) + V 2 [(V * 0)(x, t)] + V • fj(x, t) + 1/2 V 2 [V p (x - x p (t)], (17) 

where x p (t) = v p t is the position of the probe. Conversely, the force exerted by the bath particles on the 
probe is 

F(t) = -Po 1/2 V [(V p * 0)(x p (t), t)} . (18) 

We are interested in the average stationnary solution for the field in the reference frame of the particle, 

</>*(x) = ($(x + Vpt, t)); (19) 


it satisfies 

T( Vp • V) 2 tf>* (x) - Vp • Vc/>*(x) = TV 2 f(x) + V 2 (V * (x)) + p“ 1/2 V 2 Vp(x). 

In Fourier space, this equation reads 

[k 2 {T + V(/c)) - iv p • k - r( v p ■ k) 2 ] (p*(k) = -p^ 1/2 /rV p (/c), 


leading to the solution 

4>*(k) 

Finally, the density variation <5p*(x) = 

Sp*(k) 

Po 


= _ -P 0 1 / 2 k 2 Vp(k) _ 

k 2 (T + V(k)) - /Vp • k - r(i/ p • k) 2 

Pg 1 ^ 2 0*(x) is, in Fourier space, 

__ ~k 2 V p (k) _ 

k 2 (T + p 0 V(k)) - iv p • k - r( v p ■ k) 2 


( 20 ) 

( 21 ) 

( 22 ) 


(23) 
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4.2 Drag coefficient 

The drag force can be expressed with the density variation in Fourier space as 


F = -Po V2 


= 'Po 1 


rJ k 

/*V P (W(*)(2^. 

kk 2 V p (k) 2 


dk 


= “Po 1 


k 2 (T + V{k)) - tv p • k - t(v p ■ k) 2 {2ix) d ' 
f k(v p ■ k)k 2 V p (k) 2 dk 


W 


[k 2 (T + V(k)) - r(v p ■ k) 2 f + {vp ■ k) 2 ( 2?r ) 

Decomposing the wavevector as k = (/cy, /c_l) according to the velocity v p allows to write the force as 


(24) 

(25) 

(26) 


-PoS 


k}k 2 Vp{kf 


dk 


[k 2 (T + V(k)) - r(v p /C||) 2 ] 2 + (v p k {{ ) 2 ( 2?r ) 
The drag coefficient A is defined by F — —Av p ; it reads 


d ‘ 


(27) 


A — po 


k 2 k 2 Vp{kf 


dk 


[k 2 {T + p 0 V(k )) - r(v p /(||) 2 ] 2 + (v p /<||) 2 ( 27r ) c/ 


(28) 


This expression is our main result. Note that this drag coefficient is induced by the bath colloids only, the 
total drag force also includes the drag induced by the solvent. 

We are interested in the stationnary bath density around the probe (Eq. 23) and the drag coefficient 
(Eq. 28) in two cases: 


(i) the probe is much larger than the bath particles, 


(ii) the probe is identical to the bath particles. 

In the case (i), the wavevectors k that contribute to the integral in Eq. 28 are of order k ~ l/a p <C 1. 
At this scale, the Fourier transform of the bath pair potential is almost constant, V(k) ~ \/(0) for k ~ l/a p . 
This allows to define the “collective” diffusion coefficient of the bath [13], 

Dbath = T + p 0 \/(0), (29) 


that sets the relaxation time of the density field p(x, t). The Peclet number is defined as 

3pV p 


Pe = 


D 


bath 


it compares the velocity of the probe to the relaxation of the density field [13]. 
The density around the probe can then be written 

5p*(k) _ _L_ -k 2 Vp(k) _ 

Po Dbath 1(2 _ j^l (|| _ rDbath^jT^I 


(30) 


(31) 


To rescale the lengths to the probe size a p , we introduce q — a p k and V p (y) — V p (a p y ) (so that V p (k) — 
a p V p (apk)). With these notations, the density around the probe reads 

£p*(q/a p) = _ -q%(q) (32) 

Po Obath q 2 - /Pe< 7 || - ^Pe 2 ^' 


5 














-20 -10 0 
Xj_/a p 


10 20 



0.1 


- 0 




Q. 


- 0.1 


x±/a p 


Figure 2: (color online) Bath density variations in the reference frame of the probe, 5p*(x ) = p*(x) — po, for 
different values of the Peclet number, when the probe is much larger than the bath particles and the bath 
particles are heavy, r/r re | = 100. The normalization constant is AT = e p /Db a th- 


where 



(33) 


is the relaxation time of the field on the lengthscale of the probe. With the appropriate normalization, this 
density depends only on the Peclet number and the inertial number r/r re |. This is also true for the drag 
coefficient, which can be written 


Po4 r qf\q 2 %{q) 2 _ dq 

^ bath ^ q 2 - ^Pe 2 g 2 ^ + Pe 2 q| ^ 7r ) 


(34) 


4.3 Numerical computation 

The dimension is set to cl — 3 and the pair potentials are Gaussian, 

(35) 

(36) 


= £ exp [ 


V(x) 

Vp(x) = e p exp I - 


2a 2 


The bath density is plotted on Fig. 2 for heavy bath particles (r/r re i = 100) at different Peclet numbers. 
The drag coefficient is plotted as a function of the Peclet number on Fig. 3, for different values of the inertial 
number. These curves exhibit shear-thickening and resemble those obtained by numerical simulations in [15] 
(see Fig. 2a). This effect is due to inertia, as it disappears at low inertial numbers, in accordance with the 
conclusions of [15]. The density profiles differ strongly from those found in [21, 18, 35, 27] without inertia, 
or in experiments [12, 13]. In the experiments of [12], one can estimate r/(a^/Do) — 10 (Do = k^T/X is 
the thermal diffusion coefficient) so that inertia should matter; however, the bath particles are hard and the 
model discussed here may not apply. Inertia allows density waves to propagate through the bath leading to 
the cone visible for Pe = 0.1, 0.5 on Fig. 2. 

In the case (ii), the dynamics of the bath depends on the mode k and no collective diffusion coefficient 
can be defined. The bath density is plotted on Fig. 4 for poe = 1, r = 100, and different probe velocities; its 
structure is more complex than in the case (i), because the bath is dispersive at the scale of the colloids. The 
drag coefficient is plotted as a function of the probe velocity on Fig. 5 for po e = 1 and different inertial times 
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Figure 3: (color online) Friction coefficient when the probe is much larger than the bath particles as a 
function of the Peclet number for different values of the inertial number r/r re |, normalized by its value at 
small Peclet number. The dots indicate the values used for the density profiles in Fig. 2. Inset: Normalized 
friction coefficient as a function of the Peclet number rescaled by \Jr/r re \. 

r, and on Fig. 6 for r = 0.1 and different interaction strengths poe. A shear-thickening regime emerges at 
large inertial times and large interaction strengths, consistently with our finding for point-like bath particles. 

4.4 Timescales and scaling laws 

We show that the different observed behaviors can be rationalized by comparing the different timescales 
involved in our process. We focus on the case (ii) where the probe is identical to the bath particles, 
a p = a — 1. Four timescales emerge in our analysis: 

• The inertial timescale T\ = m/ A s = r. 

• The thermal diffusion timescale r t h = A s a p /T = 1, 

• The timescale associated with the density relaxation due to pair interactions, r pa ■„ = A s /(poa p e) — 
l/(poe). The denominator, p^a^e, is the energy scale seen by one bath particle. This expression differs 
slightly from the one given in [15] for the viscous damping timescale (tq in [15]), where po ~ a~ 3 and 
T"pair = ^s3 p /e. 

• The timescale associated with the motion of the probe, r p = a p /v p = 1/V P . In macrorheology, the 
characteristic timescale of the forcing is set by the shear rate 7 : T S h ear = 1 / 7 . 

At moderate density and low temperature, the suspension undergoes a glass transition [15]. We are 
interested in higher densities and temperatures, where the system still behaves as a fluid [40, 41, 42, 43]. 
In this case, thermal diffusion and pair interactions act together to relax the density; the first two terms on 
the r.h.s of Eq. 16 show that the two associated timescales r t h and r pa j r are combined in one timescale r re 1 
related to the relaxation of the density field: 

1 1 1 

— —-!-— T + Pq6. (37) 

Tel Th r pair 
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Figure 4: (color online) Bath density variations in the reference frame of the probe, p*(x)/po — 1, for different 
values of the probe velocity when the bath particles are identical to the probe. The bath particles are heavy, 
r = 100, and poe = 1. 
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Figure 5: (color online) Friction coefficient as a function of the probe velocity when the bath particles are 
identical to the probe, for different values of the inertial relaxation time r, with poe = 1. The friction 
coefficient is normalized by its value at small velocity. The dots indicate the values used for the density 
profiles in Fig. 4. Inset: Normalized friction coefficient as a function of the probe velocity rescaled by y/r. 
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Figure 6: (color online) Friction coefficient as a function of the probe velocity when the bath particles 
are identical to the probe, for different values of the interaction strength poe, with r = 0.1. The friction 
coefficient is normalized by the density of the bath. 

Without inertia, i.e. T\ = 0, only two regimes are present (see Fig. 5, r — 0.1): a Newtonian regime 
at small velocities where r p > T re |, and a shear-thinning regime at large probe velocities such that r p < 
r re |. The suspension shear thins because it does not have enough time to respond to the presence of the 
probe. This effect has already been evidenced in this framework in [27] and for a probe in more general 
environments in [39]. Shear-thinning is also observed in experiments [12, 13], numerical simulations [18, 15] 
and computations [20, 21, 24, 25, 23, 26]; it is due to the disruption of the static organization of the medium 
by the probe, which occurs when the forcing timescale becomes comparable to the relaxation timescale. 

When inertia is added, expanding the denominator of the integrand in Eq. 28 for small probe velocities 
leads to the following criterion: shear-thickening is present if 

T>T r el- (38) 

This is also the condition for the mode k = 7r/a p of the density to oscillate and waves to develop in the 
wake of the probe (see Figs. 2 and 4). This expansion gives the variation of the friction coefficient at small 
velocity as 

A(v p )-A(v p = 0) ~ f-^-l)v p 2 ; (39) 

Vp->0 \T re | J 

this scaling law has been observed in [15], but it does not hold when the maximal drag coefficient is approached 
here. It is also different from the scaling A ~ v p obtained by Bagnold [6], where the dissipation is dominated 
by the collisions between the grains. Note that the exponent is dictated by symmetry in our model: the drag 
coefficient is an even function of v p , so that the first correction is of order v p . 

A direct look at the same expression shows that the viscosity is maximal when 

t p ^ VTTei, (40) 

which corresponds to a resonance between the forcing and the excited mode. This expression does not 
match the scaling found in [15] and explained by arguments from kinetic theory, which is r p ~ ’ T \^ T ]l\ r ■ 
Flowever, Kawasaki et al. suggest that the scaling (40) can be obtained assuming “soft particles and collisional 
dissipation”. 
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As the inertia r increases, the shear-thickening becomes sharper in Figs. 3 and 5. Rescaling the velocity 
by s/r in Eq. 28 (or the Peclet number by y/r/r re \ in Eq. 34), it appears that the integral diverges for 
t p > V T i r rei when r —» oo. This divergence is shown in the inset of Figs. 3 and 5; it shows that the 
shear-thickening becomes discontinuous in this limit. 

At high velocities such that r p <C ^T\T re \, the drag coefficient induced by the bath particles decays as 
~ i/“ 4 . The total drag force thus reduces to the drag force induced by the solvent. 


4.5 Constant velocity vs. constant force 

As pointed out above, the drag coefficient depends on the mode of driving: constant force or constant 
velocity [35, 14], We show this effect in the limit of zero velocity, or zero force, and without inertia (r = 0). 
The drag coefficient (28) reduces to 


Po f V(k) 2 dk 
d J [7- + p 0 \/(/c)] 2 (27r) rf ' 


The drag coefficient at constant force has been computed in [27] (Eq. 79); at zero force, it is 


A Po f VW 2 _ c ^_ (A2) 

°" d J [T + p 0 V{k)}[ 2T + p 0 \/(/c)] (27r) rf - 

As found in [35], the drag coefficient is smaller when measured at constant force, because of the factor 2 
in front of the temperature in the denominator of Eq. 42, which accounts for the diffusion of the probe [27], 
Interestingly, in the “dilute” limit where poV(k) <C T, we recover the relation A cv . = 2A c .f. found in [35]. 


5 Conclusion 

We considered a suspension of colloids with inertia and showed that the drag force on a probe pulled at 
constant velocity can be computed analytically in the limit of a dense suspension and soft colloids. The 
rheology of the suspension can be deduced: shear-thickening is observed if the inertia is large enough for 
density waves to propagate through the solution, in agreement with [15], However, some quantitative 
differences arise between these numerical simulations and our computations in the scalings of the viscosity 
before the maximum and the position of the maximum. Numerical simulations closer to the regime that we 
studied, which are beyond the scope of this paper, would help to bridge this gap. 

In [27], we shown that the ideal limit of a very dense and soft suspension reproduces the tracer diffusion 
coefficient computed in the dilute limit [44] and exhibits force-induced diffusion, which has also been observed 
for hard particles [45, 23, 46]. Adding inertia, we showed here that it reproduces the basic features of inertia- 
induced shear-thickening reported for numerical simulations for hard particles [15]. These examples show 
that albeit the assumptions made in the computation are restrictive, this mean-field model can shed light on 
generic rheological properties of colloidal suspensions. At this stage, when applied to systems that are not in 
the theoretical range of validity, the outcome of this model cannot be considered as quantitatively accurate. 
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